El objetivo de nuestro trabajo es estudiar si existe alguna relación entre la vacuna BCG (Bacillus de Calmette y Guérin) para la tuberculosis y los datos de mortalidad de la COVID-19 en algunos países, ya que hay estudios que sugieren esta vacuna incrementa las capacidades inmunitarias de la población, hecho que se ve en el número reducido de fallecimientos por COVID-19 en ciertos países. Mediante los conjuntos de datos de BCG y de mortalidad por COVID-19 cedidos por The BCG world atlas y por BCG - COVID-19 AI Challenge de Kaggle, vamos a intentar desvelar dichas relaciones.
Los ficheros en cuestión son del tipo `csv, así que son fácilmente importables a data frames enR:
# Cargamos ambos datasets,
BCG_strain <-
read_csv("task_2-BCG_strain_per_country-1Nov2020.csv")
COVID_noformat <-
read_csv(
"task_2-COVID-19-death_cases_per_country_after_fifth_death-till_22_September_2020.csv"
)
# Intenté ver que hay dentro de los data frames, pero el print es feo así que lo escribiré a mano
# str(COVID_noformat)
# str(BCG_strain)
El contenido de las variables BCG_strain y COVID_noformat es el siguiente:
| BCG_strain | COVID_noformat |
|---|---|
| country_name | country_name |
| country_code | alpha_3_code |
| mandatory_bcg_strain_2015-2020 | date_first_death |
| mandatory_bcg_strain_2010-2015 | date_fifth_death |
| mandatory_bcg_strain_2005-2010 | deaths_per_million_10_days_after_fifth_death |
| mandatory_bcg_strain_2000-2005 | deaths_per_million_15_days_after_fifth_death |
| mandatory_bcg_strain_1990-2000 | deaths_per_million_20_days_after_fifth_death |
| mandatory_bcg_strain_1980-1990 | deaths_per_million_25_days_after_fifth_death |
| mandatory_bcg_strain_1970-1980 | deaths_per_million_30_days_after_fifth_death |
| mandatory_bcg_strain_1960-1970 | deaths_per_million_35_days_after_fifth_death |
| mandatory_bcg_strain_1950-1960 | deaths_per_million_40_days_after_fifth_death |
| vaccination_timing_unified | deaths_per_million_45_days_after_fifth_death |
| BCG Atlas: Which year was vaccination introduced? | deaths_per_million_50_days_after_fifth_death |
| Year of changes to BCG schedule | deaths_per_million_55_days_after_fifth_death |
| BCG Atlas: BCG Recommendation Type | deaths_per_million_60_days_after_fifth_death |
| BCG Atlas: Details of changes | deaths_per_million_65_days_after_fifth_death |
| BCG Atlas: Timing of 1st BCG? | deaths_per_million_70_days_after_fifth_death |
| BCG Atlas: BCG Strain | deaths_per_million_75_days_after_fifth_death |
| BCG Atlas: How long has this BCG vaccine strain been used? | deaths_per_million_80_days_after_fifth_death |
| deaths_per_million_85_days_after_fifth_death | |
| deaths_per_million_90_days_after_fifth_death | |
| deaths_per_million_95_days_after_fifth_death | |
| deaths_per_million_100_days_after_fifth_death | |
| deaths_per_million_105_days_after_fifth_death | |
| deaths_per_million_110_days_after_fifth_death | |
| deaths_per_million_115_days_after_fifth_death | |
| deaths_per_million_120_days_after_fifth_death | |
| deaths_per_million_125_days_after_fifth_death | |
| deaths_per_million_130_days_after_fifth_death | |
| deaths_per_million_135_days_after_fifth_death | |
| deaths_per_million_140_days_after_fifth_death | |
| deaths_per_million_145_days_after_fifth_death | |
| deaths_per_million_150_days_after_fifth_death | |
| stringency_index_10_days_after_fifth_death | |
| stringency_index_15_days_after_fifth_death | |
| stringency_index_20_days_after_fifth_death | |
| stringency_index_25_days_after_fifth_death | |
| stringency_index_30_days_after_fifth_death | |
| stringency_index_35_days_after_fifth_death | |
| stringency_index_40_days_after_fifth_death | |
| stringency_index_45_days_after_fifth_death | |
| stringency_index_50_days_after_fifth_death | |
| stringency_index_55_days_after_fifth_death | |
| stringency_index_60_days_after_fifth_death | |
| stringency_index_65_days_after_fifth_death | |
| stringency_index_70_days_after_fifth_death | |
| stringency_index_75_days_after_fifth_death | |
| stringency_index_80_days_after_fifth_death | |
| stringency_index_85_days_after_fifth_death | |
| stringency_index_90_days_after_fifth_death | |
| stringency_index_95_days_after_fifth_death | |
| stringency_index_100_days_after_fifth_death | |
| stringency_index_105_days_after_fifth_death | |
| stringency_index_110_days_after_fifth_death | |
| stringency_index_115_days_after_fifth_death | |
| stringency_index_120_days_after_fifth_death | |
| stringency_index_125_days_after_fifth_death | |
| stringency_index_130_days_after_fifth_death | |
| stringency_index_135_days_after_fifth_death | |
| stringency_index_140_days_after_fifth_death | |
| stringency_index_145_days_after_fifth_death | |
| stringency_index_150_days_after_fifth_death |
Una visualización preliminar de estos datos revela que son todos del tipo string y que además muchas columnas sin datos (columnas cuyo único contenido es NULL), por lo tanto llevaremos a cabo una limpieza de los mismos además de cambios de tipo de variables para que las manipulaciones posteriores sean más cómodas. Los detalles se muestran en el siguiente bloque de código:
# Limpiar datos de BCG
# Elimino columnas que sean sólo NA
BCG_strain <- BCG_strain[, apply(!is.na(BCG_strain), 2, all)]
# De momento, no me interesa qué vacunas se ponían cada año, sino si se ponían o no.
# Transformo los valores de cada año en
# 0 - No se ponía vacuna, hasta ahora None
# 1 - Sí se ponía vacuna
# NA - Este dato es desconocido, hasta ahora Unknown
BCG_strain_no_strain <- BCG_strain
# Transformo los valores de las columnas
BCG_strain_no_strain[, -1] <-
sapply(BCG_strain_no_strain[, -1], function(x) {
a <-
gsub("None", 0, x) %>% gsub("Unknown", NA, .) # Añado los 0 y los NA.
for (i in 1:length(a)) {
# Serán 1 aquellos que no sean ni 0 ni NA
if (a[i] != "0" && !is.na(a[i])) {
a[i] <- 1
}
}
return(as.integer(a)) # Cambio las columnas a integer
})
BCG_no_strain_no_NA <- na.omit(BCG_strain_no_strain)
#Versiñon más compacta del dataframe, sin datos a diferentes días o años.
#Agrupando los datos de vacunas en tres columnas:
#periods_with_vaccine - incluye el número de periodos estudiados con vacunación activa
#vaccination_2020_2015 - el único periodo con el que nos quedamos, el último
#first_vaccine_year - de los años estudiados, el primero con campaña de vacunación. En el caso de no tener vacunación, este será el último año estudiado (2020)
#last_vaccine_year - de los años estudiados, el último con campaña de vacunación. En el caso de no tener vacunación, este será el primer año estudiado (1950)
#Creamos el nuevo dataframe simplificado
BCG_no_strain_simple <- data.frame(
"country_name" = BCG_no_strain_no_NA$country_name,
"periods_with_vaccine" = BCG_no_strain_no_NA%>%
.[2:ncol(.)] %>% rowSums(), #sumamos los periodos con vacuna
"vaccination_2020_2015" = BCG_no_strain_no_NA$`mandatory_bcg_strain_2015-2020`)
#Añadimos el último año con vacunación
BCG_no_strain_simple$last_vaccine_year <- names(BCG_no_strain_no_NA[2:ncol(BCG_no_strain_no_NA)])[max.col(BCG_no_strain_no_NA[2:ncol(BCG_no_strain_no_NA)] != 0, 'first')] %>%
substring(nchar(.)-3, nchar(.))%>%
as.numeric()
#Añadimos el primer año con vacunación
BCG_no_strain_simple$first_vaccine_year <- names(BCG_no_strain_no_NA[2:ncol(BCG_no_strain_no_NA)])[max.col(BCG_no_strain_no_NA[2:ncol(BCG_no_strain_no_NA)] != 0, 'last')] %>%
substring(nchar(.)-8, nchar(.)-5)%>%
as.numeric()
#El próximo código es necesario para que los países sin campaña de vacunación no obtengan los mejores valores. El código utilizado para obtener el último o primer año de vacunación les favorece, ya que obtiene el primer o el último índice de aquellos valores distintos de 0. Como en su caso no hay ningún valor distinto a 0, este sería simplemente el primero o el último. Para arreglar esto, establezco manualmente que tengan el último año de vacunación más bajo posible y el primer año de vacunación más alto posible.
BCG_no_strain_simple[BCG_no_strain_simple$periods_with_vaccine == 0,]$last_vaccine_year = 1950
BCG_no_strain_simple[BCG_no_strain_simple$periods_with_vaccine == 0,]$first_vaccine_year = 2020
####################################################################################
# Limpiar datos de COVID
# Elimino columnas que sean sólo NA
COVID_noNA <- COVID_noformat[, apply(!is.na(COVID_noformat), 2, all)]
# En este caso, para variar, los valores vacíos están denotados como NULL,
# cambio esto a NA
COVID_Na <- sapply(COVID_noNA, function(x)
gsub("NULL", NA, x))
# El resulatado de la función anterior es una string. Lo convierto a dataframe.
COVID_Na_df <- as.data.frame(COVID_Na)
# Modifico las fechas para que se almacenen como Date
COVID_Na_df[, c("date_fifth_death")] <-
as.Date(COVID_Na_df[, c("date_fifth_death")], "%d/%m/%y")
COVID_Na_df[, c("date_first_death")] <-
as.Date(COVID_Na_df[, c("date_first_death")], "%d/%m/%y")
# Modifico las muertes para que se almacenen como floats.
COVID_Na_df[, -c(1, 2, 3, 4)] <-
sapply(COVID_Na_df[, -c(1, 2, 3, 4)], as.numeric)
COVID_Na_df2 <- data.frame("country_name" = COVID_Na_df$country_name,
"dpm_100d" = COVID_Na_df$deaths_per_million_100_days_after_fifth_death,
"si_100d" = COVID_Na_df$stringency_index_100_days_after_fifth_death)
################################################################################################
# Junto ambos dataframes en uno sólo.
COVID_BGC <-
inner_join(BCG_no_strain_simple, COVID_Na_df2, by = "country_name")
Nuestra tabla resultante es la siguiente:
| country_name | periods_with_vaccine | vaccination_2020_2015 | last_vaccine_year | first_vaccine_year | dpm_100d | si_100d |
|---|---|---|---|---|---|---|
| Afghanistan | 6 | 1 | 2020 | 1980 | 26.00 | 78.70 |
| Angola | 6 | 1 | 2020 | 1980 | 4.38 | NA |
| Argentina | 3 | 0 | 2010 | 1990 | 30.64 | 92.59 |
| Armenia | 4 | 1 | 2020 | 2000 | 190.67 | NA |
| Australia | 9 | 1 | 2020 | 1950 | 4.00 | 52.31 |
| Bangladesh | 6 | 1 | 2020 | 1980 | 11.95 | 74.54 |
Mediante esta tabla llevaremos a cabo nuestros análisis. A continuación mostramos la estructura de la misma:
str(COVID_BGC)
## 'data.frame': 65 obs. of 7 variables:
## $ country_name : chr "Afghanistan" "Angola" "Argentina" "Armenia" ...
## $ periods_with_vaccine : num 6 6 3 4 9 6 6 9 8 9 ...
## $ vaccination_2020_2015: int 1 1 0 1 1 1 1 1 1 1 ...
## $ last_vaccine_year : num 2020 2020 2010 2020 2020 2020 2020 2020 2020 2020 ...
## $ first_vaccine_year : num 1980 1980 1990 2000 1950 1980 1980 1950 1960 1950 ...
## $ dpm_100d : num 26 4.38 30.64 190.67 4 ...
## $ si_100d : num 78.7 NA 92.6 NA 52.3 ...
a) ¿Cuál es la media para la variable “deaths per million 150 days after fifth death” (representado por `dpm_100d) para el conjunto de países que la han medido? ¿Y para la variable “stringency index 150 days after fifth death”?
m1 <- mean(na.omit(COVID_BGC$dpm_100d))
m2 <- mean(na.omit(COVID_BGC$si_100d))
cat(sprintf("La media para dpm_100d es %.2f\n", m1))
## La media para dpm_100d es 110.44
cat(sprintf("La media para si_100d es %.2f\n", m2))
## La media para si_100d es 57.84
b) ¿Cuáles son los países cuyo valor para “deaths per million 100 days after fifth death” (en caso de estar presente) se encuentra por debajo de la media? ¿Y para la variable “stringency index 100 days after fifth death”?
cat(
"Los países cuyo valor de dpm_100d es menor que la media del dataset son:\n\n",
paste(subset(COVID_BGC, dpm_100d < m1)$country_name, collapse = ', '),
"\n\n"
)
## Los países cuyo valor de dpm_100d es menor que la media del dataset son:
##
## Afghanistan, Angola, Argentina, Australia, Bangladesh, Bosnia and Herzegovina, Bulgaria, Colombia, Czech Republic, Denmark, El Salvador, Estonia, Ethiopia, Finland, Greece, Guam, Hungary, India, Indonesia, Japan, Kazakhstan, Kenya, Kuwait, Latvia, Malaysia, Malta, Nepal, Nigeria, Pakistan, Philippines, Poland, Saudi Arabia, Senegal, Sierra Leone, South Africa, Sudan, Taiwan, Tanzania, Thailand, Tunisia, Ukraine, Uruguay
cat(
"Los países cuyo valor de si_100d es menor que la media del dataset son:\n\n",
paste(subset(COVID_BGC, si_100d < m2)$country_name, collapse = ', ')
)
## Los países cuyo valor de si_100d es menor que la media del dataset son:
##
## Australia, Bosnia and Herzegovina, Bulgaria, Czech Republic, Estonia, Finland, Greece, Hungary, Ireland, Italy, Japan, Latvia, Malaysia, Poland, Senegal, Sierra Leone, Spain, Sweden, Switzerland, Taiwan, Tanzania, Thailand, Tunisia, Ukraine, Uruguay
c) ¿Cuáles son los países que cumplen ambas condiciones del apartado anterior?
cat(
"Los países cuyos valores de dpm_100d y de si_100d son menores que la media del dataset son:\n\n",
paste(subset(COVID_BGC, (dpm_100d < m1) &
(si_100d < m2))$country_name, collapse = ', ')
)
## Los países cuyos valores de dpm_100d y de si_100d son menores que la media del dataset son:
##
## Australia, Bosnia and Herzegovina, Bulgaria, Czech Republic, Estonia, Finland, Greece, Hungary, Japan, Latvia, Malaysia, Poland, Senegal, Sierra Leone, Taiwan, Tanzania, Thailand, Tunisia, Ukraine, Uruguay
d) ¿Cuáles son los países que han tenido campaña de vacunación de la vacuna BCG más reciente y que su la media de mortalidad a los 150 días es menor que la media?
cat(
"Los países cuyos valores de dpm_100d son menores que la media del dataset y que además han tenido una reciente campaña de vacunación de la vacuna BCG son:\n\n",
paste(subset(
COVID_BGC, (dpm_100d < m1) &
(vaccination_2020_2015 == 1)
)$country_name, collapse = ', ')
)
## Los países cuyos valores de dpm_100d son menores que la media del dataset y que además han tenido una reciente campaña de vacunación de la vacuna BCG son:
##
## Afghanistan, Angola, Australia, Bangladesh, Bosnia and Herzegovina, Bulgaria, Colombia, El Salvador, Estonia, Ethiopia, Hungary, India, Indonesia, Japan, Kazakhstan, Kenya, Kuwait, Latvia, Malaysia, Malta, Nepal, Nigeria, Pakistan, Philippines, Poland, Saudi Arabia, Senegal, Sierra Leone, South Africa, Sudan, Taiwan, Tanzania, Thailand, Tunisia, Ukraine
a) Resumen estadístico de las variables
# Obvimente hay variables en las que no tiene sentido hacer resumen estadístico, como el alpha_3_code, las strains... Pero por ahora lo voy a dejar
summary(COVID_BGC)
## country_name periods_with_vaccine vaccination_2020_2015
## Length:65 Min. :0.000 Min. :0.0000
## Class :character 1st Qu.:5.000 1st Qu.:1.0000
## Mode :character Median :6.000 Median :1.0000
## Mean :6.277 Mean :0.7538
## 3rd Qu.:8.000 3rd Qu.:1.0000
## Max. :9.000 Max. :1.0000
##
## last_vaccine_year first_vaccine_year dpm_100d si_100d
## Min. :1950 Min. :1950 Min. : 0.294 Min. :19.44
## 1st Qu.:2020 1st Qu.:1950 1st Qu.: 9.524 1st Qu.:38.89
## Median :2020 Median :1960 Median : 30.134 Median :59.72
## Mean :2012 Mean :1969 Mean :110.444 Mean :57.84
## 3rd Qu.:2020 3rd Qu.:1980 3rd Qu.:116.554 3rd Qu.:74.07
## Max. :2020 Max. :2020 Max. :584.686 Max. :92.59
## NA's :9 NA's :12
# Igual habría que dibujar histogramas de tasas de muerte, y de vacunas puestas. No sé la verdad.
b) Distribución de las variables
Observemos cómo se distribuyen los datos de “deaths per million 100 days after fifth death” y de “stringency index 100 days after fifth death”
ggplot(COVID_BGC, aes(x= dpm_100d, fill = "b"))+geom_density()+
theme_bw()+theme(legend.position = 0, panel.grid = element_blank())+
ggtitle("Distribución de muertes por millón a los 100 días de la quinta muerte")
ggplot(COVID_BGC, aes(x= si_100d, fill = "b"))+geom_density()+
theme_bw()+theme(legend.position = 0, panel.grid = element_blank())+
ggtitle("Distribución del stringency index por millón a los 100 días de la quinta muerte")
Mediante el gráfico de densidad podemos ver la distribución de dos variables continuas como son el número de muertes por millón y el stringency index. Podemos comprobar la distribución de una variable binomial como es la presencia o no de campaña de vacunación en 9 periodos de 5 o 10 años.
ggplot(COVID_BGC, aes(x= vaccination_2020_2015, fill = "b"))+geom_bar()+
theme_bw()+theme(legend.position = 0, panel.grid = element_blank())+
ggtitle("Presencia de una campaña de vacunación entre 2015 y 2020")
ggplot(COVID_BGC, aes(x= periods_with_vaccine, fill = "b"))+geom_bar()+
theme_bw()+theme(legend.position = 0, panel.grid = element_blank())+
ggtitle("Número de décadas o lustros con campaña de vacunación activa")
cormat <-
cor(COVID_BGC[2:ncol(COVID_BGC)]%>%na.omit())
cormat2 <- cormat
cormat2[upper.tri(cormat2)] <-
NA #Para visualizar solamente una vez las correlaciones
cormat2 <- melt(round(cormat2, 2)) #Formato para poder usar ggplot
ggplot(cormat2, aes(x = Var1, y = Var2, fill = value)) + geom_tile() + scale_fill_continuous(type = "viridis")
fig <-
plot_ly(
x = colnames(cormat),
y = colnames(cormat),
z = cormat,
type = "heatmap"
)
fig
ggplot(COVID_BGC,
aes(x = dpm_100d, y = vaccination_2020_2015, label = country_name)) +
geom_jitter(position = position_jitter(seed = 1)) +
geom_label_repel(size = 2, position = position_jitter(seed = 1)) +
xlim(c(-100, 800))
ggplot(COVID_BGC,
aes(x = dpm_100d, y = last_vaccine_year, label = country_name)) +
geom_jitter(position = position_jitter(seed = 1)) +
geom_label_repel(size = 2, position = position_jitter(seed = 1)) +
xlim(c(-100, 800))
ggplot(COVID_BGC,
aes(x = dpm_100d, y = first_vaccine_year, label = country_name)) +
geom_jitter(position = position_jitter(seed = 1)) +
geom_label_repel(size = 2, position = position_jitter(seed = 1)) +
xlim(c(-100, 800))
model <- lm(dpm_100d ~ vaccination_2020_2015 + first_vaccine_year + last_vaccine_year, COVID_BGC)
summary(model)
##
## Call:
## lm(formula = dpm_100d ~ vaccination_2020_2015 + first_vaccine_year +
## last_vaccine_year, data = COVID_BGC)
##
## Residuals:
## Min 1Q Median 3Q Max
## -319.96 -66.61 -17.95 12.25 385.72
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16054.899 4441.476 3.615 0.000679 ***
## vaccination_2020_2015 -44.886 61.859 -0.726 0.471328
## first_vaccine_year -2.499 1.024 -2.441 0.018081 *
## last_vaccine_year -5.465 1.597 -3.421 0.001221 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 127.5 on 52 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.4378, Adjusted R-squared: 0.4054
## F-statistic: 13.5 on 3 and 52 DF, p-value: 1.241e-06